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Probabilistic  Sensitivity  Analysis  of  Fretting  Fatigue 


Patrick  J.  Golden 

Air  Force  Research  Laboratory,  Wright-Patterson  AFB,  OH  45433 

Harry  R.  Millwater'  and  Xiaobin  Yang1 
University  of  Texas  at  San  Antonio,  San  Antonio,  TX  78249 


A  probabilistic  analysis  of  the  fatigue  life  of  specimens  subject  to  fretting  fatigue  was 
carried  out.  A  mechanics  based  fretting  life  analysis  was  applied  that  accounted  for  the  local 
stress  gradient  at  the  edge  of  contact.  The  random  variables  in  the  analysis  included  the 
initial  crack  size,  coefficient  of  friction,  crack  growth  rate  law,  and  the  contact  pad  profile. 
The  variation  in  pad  profiles  was  determined  through  statistical  analysis  of  seventy -seven 
machined  pads.  A  probabilistic  fatigue  analysis  was  applied  using  Monte  Carlo  sampling  to 
determine  the  statistics  (mean  and  standard  deviation)  of  the  fatigue  life  prediction  and  the 
probabilistic  sensitivities  (partial  derivatives  of  the  fatigue  statistics  with  respect  to  the  input 
probability  density  function  parameters). 


Nomenclature 


a 

= 

slope  of  shear  force  versus  normal  force  in  partial  slip 

M 

= 

coefficient  of  friction 

Gy 

= 

components  of  stress 

a 

= 

crack  depth 

c 

= 

half  surface  crack  length 

F 

= 

applied  load  in  the  fretting  test 

h(x) 

= 

contact  gap  function 

K 

= 

stress  intensity  factor 

M 

= 

contact  moment 

R 

= 

load  ratio 

P(x) 

= 

contact  pressure  traction 

P 

= 

normal  contact  force 

q(x) 

= 

contact  shear  traction 

Q 

= 

shear  contact  force 

I.  Introduction 

Fretting  is  a  problem  in  many  aerospace  applications  including  the  blade  to  disk  attachment  in  turbine  engines. 
Two  fretting  modes  can  contribute  to  damage  in  fretting:  gross  slip  when  the  two  surfaces  slide  resulting  in  wear, 
and  partial  slip  when  the  two  surfaces  are  nominally  stuck  together  except  for  a  small  slip  zone  at  the  edge  contact. 
The  surface  damage  and  wear  caused  by  fretting  is  a  costly  maintenance  burden  and  when  combined  with  the  very 
high  local  contact  stresses  due  to  fretting,  it  can  result  in  disk  or  blade  cracking  and  the  potential  for  catastrophic 
failure.  A  probabilistic  analysis  is  often  applied  to  problems  to  help  determine  the  effect  of  variability  of  the  model 
input  parameters  on  the  model  outputs.  Prior  work  on  statistical  or  probabilistic  analysis  in  fretting  includes 
modeling  of  the  variability  in  the  contact  surface  profile  by  Kumari  and  Farris  [1],  Here,  the  measured  profiles  of 
the  indentors  in  fretting  fatigue  tests  were  measured  and  statistically  described  and  carried  through  the  life  prediction 
models  to  estimate  the  expected  variability  in  stress  and  life.  Other  probabilistic  fretting  analyses  include  work  on 
fretting  fatigue  of  riveted  lap  joints  [2]  and  on  rolling  contact  [3],  The  objective  of  the  current  work  is  to  develop 
and  demonstrate  a  probabilistic  fretting  fatigue  lifing  approach  for  a  dovetail  fretting  experimental  configuration  that 
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includes  a  broad  range  of  input  random  variables,  and  to  apply  efficient  sensitivity  methods  to  determine  the  relative 
importance  of  the  input  variables.  A  previously  demonstrated  stress  and  life  prediction  analysis  is  adapted  to  this 
effort.  Input  probability  density  functions  (PDF’s)  were  developed  using  available  laboratory  data.  A  modular 
probabilistic  sensitivity  code  was  developed  and  applied  to  this  analysis. 


Figure  1.  Schematic  of  the  dovetail  experiment  and  the 
analysis. 


II.  Deterministic  Fretting  Analysis 

A  series  of  fretting  experiments  was  previously  conducted  [4]  to  improve  understanding  of  fretting  behavior  in 
Ti-6A1-4V  and  to  test  life  prediction  models.  The  geometry  of  the  fretting  samples  was  a  dovetail  shaped  specimen 
that  was  designed  to  represent  the  attachment 

between  a  turbine  engine  blade  and  disk.  The  tests  Experiment  Analysis 

were  conducted  at  room  temperature,  which  is 
consistent  with  the  operating  conditions  of  a  fan 
disk.  The  contact  interface  was  bare  Ti-6A1-4V  on 
bare  Ti-6A1-4V,  but  several  coatings  and  residual 
stress  surface  treatments  applied  to  the  contact 
interface  were  also  tested.  The  analysis  in  the 
current  work  was  limited  to  the  bare  Ti-6A1-4V 
tests.  A  schematic  of  the  fretting  fatigue  test  rig 
modeled  in  this  analysis  is  shown  in  Figure  1 .  The 
experimental  setup  consists  of  the  dovetail 
specimen  (one-half  of  the  specimen  is  shown  in 
schematic),  two  fretting  pads,  and  a  steel  fixture. 

The  fretting  pads  are  held  in  the  steel  fixture  at  a 

45°  flank  angle,  and  the  dovetail  specimen  is  pulled  with  a  cyclic  load,  F ,  into  the  fretting  pads.  Both  the  dovetail 
specimen  and  pads  were  machined  from  Ti-6A1-4V  with  a  thickness  of  7.62  mm,  modulus  E  =  116  GPa,  and 
Poisson’s  ratio  v=  0.31.  The  nominal  fretting  pad  geometry  was  a  3  mm  flat  with  3  mm  blending  radii.  The 
normal,  P,  and  shear,  Q,  contact  forces  were  measured  indirectly  from  strain  gage  instrumentation  of  the 
experiments.  More  details  of  the  experimental 
setup  and  instrumentation  can  be  found  in  Golden 
and  Nicholas  [4].  10  tests  were  run  at  values  of 
Fmax  ranging  from  18  kN  to  30  kN  at  a  load  ratio  R 
=  0.1. 

A  fracture  mechanics  based  fretting  life 
prediction  analysis  for  dovetail  specimens  was 
developed  and  demonstrated  using  these  and  other 
experimental  test  data  and  is  described  in  Golden 
and  Calcaterra  [5].  The  objective  of  this  prior 
work  was  to  demonstrate  lifing  methods  that  could 
be  applied  to  more  complex  3-D  structures  found 
in  turbine  engines,  however,  the  current  problem 
can  be  simplified  to  a  2-D  problem.  The  analysis 
was  broken  into  two  parts:  a  finite  element  method 
(FEM)  analysis  and  a  2-D  numerical  contact  stress 
analysis.  A  finite  element  analysis  is  needed  in  this 
problem  to  determine  two  sets  of  quantities.  The 
first  quantity,  described  here,  relates  to  the 
behavior  of  the  contact  force  history.  The  second, 
which  can  be  obtained  from  the  same  FEM  model, 
is  the  bulk  stress  as  shown  in  Figure  1  and  is 
described  below.  The  dovetail  geometry  was 
modeled  using  a  nonlinear  contact  FEM  model. 

This  analysis  yields  the  contact  force  history  for  P 
and  Q  and  an  example  is  plotted  in  Figure  2.  This 

example  shows  a  single  loading  and  unloading  starting  from  no  load.  As  the  load  is  applied,  the  contact  is  initially 
in  sliding  (gross  slip)  and  the  contact  forces  follow  the  dashed  line  defined  by  the  equation 


P  (N/mm) 

Figure  2.  A  typical  plot  of  Q  versus  P  showing  the  partial 
slip  slope  a  and  the  boundaries  defined  by  /j. 
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Q  =  nP 


(i) 


where  //  is  the  coefficient  of  friction.  Upon  load  reversal,  the  contact  load  path  changes  slope  to  the  partial  slip 
slope,  a ,  which  is  a  characteristic  of  the  component  geometry  and  compliance.  The  contact  forces  will  follow  the 
slope  a  for  both  increasing  and  decreasing  loads  until  the  dashed  lines  defining  friction  are  reached.  Therefore,  the 
two  key  inputs  needed  for  prediction  of  the  contact  force  history  for  a  given  applied  load  history  are  //  and  a.  The 
coefficient  of  friction  is  a  property  that  can  only  be  measured,  however,  a  can  be  predicted  from  the  FEM  analysis 
or  measured  by  experiment.  Once  these  2  parameters  are  known,  the  contact  forces  can  be  determined  for  a  given 
applied  force  without  additional  FEM  analysis  through  a  procedure  developed  by  Gean  and  Farris  [6], 

Once  the  contact  forces  were  known,  the 
contact  stresses  were  calculated  using  a  2-D 
numerical  contact  stress  analysis.  Figure  1  shows  a 
schematic  of  the  equivalent  contact  geometry, 
defined  by  the  gap  function,  h(x),  used  in  this 
analysis.  In  the  experiments  described  here,  h(x) 
was  simply  the  profile  of  the  fretting  pads  that  are 
pressed  into  contact  with  the  flat  dovetail 
specimens.  Rather  than  use  the  prescribed  profiles 
in  the  analysis,  the  as-machined  profiles  were 
obtained  using  a  contacting  profilometer  to 
quantify  the  actual  profiles.  The  analytical  tool, 

CAPRI  (Contact  Analysis  for  Profiles  of  Random 
Indenters)  described  in  McVeigh  et  al.  [7],  was 
developed  to  solve  the  singular  integral  equation 
that  defines  this  contact  problem  using  any 
reasonable  function,  h(x),  as  measured  from  the 
profilometry.  The  output  of  CAPRI  is  the  pressure, 
p(x),  and  shear,  q(x),  tractions  due  to  the  applied 
contact  forces.  An  example  of  these  tractions  is  shown  in  Figure  3.  Next,  CAPRI  was  used  to  determine  the 
subsurface  stress  field  under  the  contact  due  to  p(x)  and  q(x).  Note  the  peaks  located  near  the  edges  of  contact  in 
Figure  3.  These  peaks  in  the  pressure,  shear,  and  the  subsurface  stresses  were  driven  by  the  geometry  near  the  edges 
of  contact  and  typically  result  in  crack  initiation  and  growth  from  the  edge  of  contact. 

The  last  step  in  the  analysis  is  the  fracture  mechanics  based  fatigue  crack  growth  prediction.  Here,  the  total 
stresses  along  the  expected  crack  growth  path  are  needed  to  calculate  the  stress  intensity  factor  range,  A K.  To  get 
the  total  stress  solution  both  the  contact  stresses  from  CAPRI  and  the  bulk  stress  as  shown  in  Figure  1  were 
superposed.  The  bulk  stress  was  the  stress  distribution  induced  in  the  component  due  to  the  remotely  applied 
loading  and  the  geometry  of  the  component.  A  procedure  to  obtain  the  bulk  stress  distribution  for  a  dovetail 
geometry  was  described  by  Golden  and  Calcaterra  [5]  which  required  both  the  FEM  model  and  CAPRI  results. 
Once  the  full  subsurface  stress  distribution  was  known,  the  mode  I  stress  gradient  could  be  extracted  along  the 
expected  crack  path  at  the  edge  of  contact.  Since  the  stress  gradient  was  nonlinear,  weight  function  stress  intensity 
factor  solutions  were  applied  [8],  The  crack  propagation  life  was  then  integrated  using  a  crack  growth  model  of  the 
form 


£  =  nm  (2) 

where  da/dN  is  the  crack  growth  rate.  Integration  of  life  was  performed  using  the  Euler  method  starting  from  a 
small  initial  fretting  crack  typically  25  mm  in  depth,  until  fracture.  Crack  growth  properties  from  previous  testing 
on  the  same  Ti-6A1-4V  material  used  in  the  US  Air  Force  High  Cycle  Fatigue  program  [9]  were  used. 

III.  Statistical  Analysis 

The  random  variables  for  the  probabilistic  fretting  analysis  performed  in  this  study  were  chosen  from  the  key 
input  variables  of  the  deterministic  analysis  described  above.  These  key  input  variables  included  the  initial  crack 
size,  the  coefficient  of  friction,  the  slope  a  that  defines  the  contact  forces  during  partial  slip,  the  crack  growth  law 
parameters,  and  the  parameters  defining  the  shape  of  the  pad  profile.  PDF’s  were  developed  for  each  of  these 
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random  variables  from  measurements  and  test  data  available  both  from  the  dovetail  fretting  experiments  and  from 
other  fatigue  and  fatigue  crack  growth  experiments.  The  details  of  the  development  of  these  PDF’s  are  described 
below.  A  summary  of  the  resulting  PDF’s  are  shown  in  Table  1. 

A.  Initial  crack  size 

The  deterministic  fretting  fatigue  life  prediction  analysis  described  above  is  based  on  fracture  mechanics  based 
crack  propagation.  To  perform  this  analysis  an  initial  crack  size  is  required.  Previous  work  [8]  has  shown  that  the 
initial  fretting  cracks  tend  to  be  shallow,  low  aspect  ratio  {ale)  surface  cracks  that  often  appear  to  be  several  shallow 
micro-cracks  that  have  coalesced.  Here,  a  is  the  surface  crack  depth  and  c  is  the  half  surface  length.  Investigations 
of  naturally  initiated  cracks  in  smooth  bar  fatigue  specimens  in  Ti-6A1-4V  with  the  same  microstructure  [10]  have 
revealed  that  the  cracks  appear  to  initiate  almost  exclusively  at  primary  alpha  grains  that  intersect  the  surface  of  the 
specimen.  These  naturally  initiated  cracks  in  the  primary  alpha  grains  form  readily  identifiable  facets  on  the 
specimen  fracture  surface.  Measurements  of  these  initial  cracks  were  made  in  a  prior  study  [10]  on  fatigue 
variability  of  Ti-6A1-4V,  where  repeated  tests  were  conducted  at  a  few  constant  amplitude  stress  levels.  The 
resulting  distribution  of  primary  alpha  grain  facets  was  used  here  to  represent  the  expected  variability  in  initial 
fretting  crack  depth.  Since  it  has  been  shown  that  typically  multiple  surface  cracks  form  and  coalesce  in  fretting,  a 
constant,  low  aspect  ratio  of  ale  =  0.2  was  used.  The  mean  and  standard  deviation  of  the  measured  naturally 
initiated  crack  sizes  is  listed  in  Table  1.  It  was  determined  that  a  lognormal  distribution  was  the  best  fit  to  the  data. 

B.  Coefficient  of  friction  and  the  partial  slip  slope 

As  described  above  and  depicted  in  Figure  2,  the  coefficient  of  friction  and  the  partial  slip  slope  of  Q  versus  P 
were  measured  during  each  test.  The  accuracy  of  these  measurements,  however,  was  limited  to  the  accuracy  of  the 
strain  gage  instrumentation  of  the  tests  and  the  finite  element  analysis  that  was  used  to  convert  fixture  strain  to 
applied  load.  Additionally,  variability  in  friction  was  expected  from  test  to  test.  Both  the  measurement  uncertainty 
and  inherent  variability  in  the  tests  needed  to  be  captured  in  the  probabilistic  model.  To  achieve  this,  data  was 
collected  from  twenty-three  fretting  tests  all  of  which  were  conducted  under  similar  loading  conditions.  All  tests 
were  bare  Ti-6A1-4V  on  Ti-6A1-4V  contact  (no  coatings)  and  all  were  loaded  at  R  =  0. 1 .  The  measured  values  of 
friction  coefficient  and  partial  slip  slope,  a,  were  found  to  be  correlated  since.  The  uncertainty  in  the  expected 
values  of  friction  and  partial  slip  slope  along  with  their  correlation  was  modeled  using  a  multivariate  normal 
distribution.  The  results  are  tabulated  in  Table  1. 

C.  Fatigue  crack  growth  parameters 

The  crack  growth  rate  curves  used  for  crack  propagation  predictions  were  based  on  fatigue  crack  growth  tests 
previously  conducted  [9],  Four  tests  from  two  different  labs  were  fit  to  a  bilinear  Paris  crack  growth  model  for  R  = 
0.1  and  -1,  where,  da/dN  is  the  crack  growth  rate  and  AK  is  the  stress  intensity  factor  range.  The  Walker  model  [11] 
was  used  to  collapse  the  data  from  the  different  R  values  using  an  equation  of  the  form 

da/dN  =c\AK(\-R),m-])\ 
da/dN  =  C2 [AK(1  -  ]"2 

where  b  is  the  intersection  point  for  regions  1  and  2,  defined  as 

^  _  l°g  io  Q  ~  l°g  io  C2 

n2  —  nt 

The  values  of  the  Walker  exponent  m  were  determined  in  the  previous  program  [9]  and  different  values  were 
used  for  R  <  0  and  R  >  0.  The  statistics  of  the  curve  fit  were  modeled  using  a  multivariate  normal.  The  coefficient 
C,  and  the  exponent  n,  were  found  to  be  correlated  for  each  segment  of  the  bilinear  crack  growth  model.  The  mean, 
standard  deviation,  and  correlation  coefficients  are  summarized  in  Table  1.  The  values  of  the  model  parameters 
resulted  units  of  MPaVm  for  A K,  and  m/cycle  for  da/dN. 


A K  <b 
AK  >  b 


(3) 


(4) 
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D.  Pad  profile 

The  nominal  geometry  for  the  contact  pad  is 
shown  in  Figure  4.  The  prescribed  profile  is  the 
geometry  that  was  requested  on  the  part  drawings 
provided  to  the  machine  shop.  The  central  segment 
of  the  prescribed  profile  is  flat  with  a  length  of  3.00 
mm.  The  edges  of  the  profile  have  constant  radii 
of  3.0  mm.  Due  to  variations  due  to  machining, 
however,  the  as-machined  profiles  differ  from  the 
prescribed  profile  in  the  part  drawing.  The  contact 
pad  profiles  were  measured  using  a  contacting 
profilometer  prior  to  being  tested.  The  contact 
profilometer  had  a  vertical  resolution  of 
approximately  10  nm.  Two  examples  of  measured 
profiles  are  plotted  in  Figure  2.  Note  that  the  scale 
on  the  y-axis  is  significantly  magnified  to  highlight 
variations  in  the  measured  profiles.  It  is  clear  from  Figure  4.  Prescribed  and  measured  machined  pad 
these  measured  profiles  that  the  central  flat  section  profiles, 
of  the  pad  is  not  truly  flat  and  that  the  radii  at  the 

edges  do  not  always  match  the  prescribed  radii  either  in  position  or  sharpness.  Additionally,  there  is  measurable 
variability  in  the  profile  geometry  from  pad  to  pad.  It  was  desired  to  capture  this  variability  through  a  statistical 
analysis  so  it  could  be  used  in  the  subsequent  probabilistic  analysis.  To  achieve  that  objective,  seventy-seven 
contact  pads  were  measured  using  the  contact  profilometer. 

In  order  to  model  the  contact  profile  variability  in  the  probabilistic  analysis,  a  mathematical  description  of  the 
contact  profile  was  needed  that  could  capture  the  important  features  of  the  profile.  Key  features  were  believed  to  be 
the  actual  radii  of  the  edges  of  the  pad  as  these  can  significantly  affect  the  stress  at  the  edge  of  contact.  Also 
important  was  the  flatness  of  the  central  portion  of  the  pad.  A  rounded  pad  versus  a  very  flat  pad  could  also 
significantly  affect  the  contact  stresses.  A  piecewise  curve  with  twenty-one  parameters  was  defined  to  represent  the 
measured  contact  profiles  and  was  written  as 


x  (|im) 


y  = 


k  x  +  b 

y°  +J(RL)~  -(x~xoJ 


kRx+bR 


x<ax 
ax  <x  <a2 
a2<  x<a3 
a3<  x<  aA 
x  >  aA 


(5) 


where  x  and  y  represent  a  Cartesian  coordinate  system  with  x  along  the  profile  and  y  the  height  of  the  profile. 
The  central  segment  of  the  profile  ( a2<x<a3 )  was  represented  by  a  sixth  order  polynomial  to  allow  the  nominally  flat 
section  to  have  curvature.  Just  outside  the  central  section  were  two  circular  arcs  to  represent  the  radii  at  the  edges  of 
contact.  Just  beyond  the  radii  were  flat  taper  sections.  Through  various  constraints  the  number  of  parameters  was 
reduced  from  twenty-one  to  thirteen.  The  constants  k  ,  k  ,  b  ,  b  ,  y0  ,  x0  ,  x0  ,  R  ,  R  ,  ch  c2,  ct,  c4,  c5,  and  c6  were 
determined  from  nonlinear  regression  of  each  of  the  measured  seventy-seven  profiles.  These  constants  were 
referred  to  as  y,  where  i  ranged  from  0  to  12.  The  remaining  8  constants  in  Equation  1  were  functions  of  y  as 
defined  by  the  constraints.  Scatter  plots  of  these  13  regression  parameters  are  shown  in  Figure  5,  which  shows  that 
many  of  the  parameters  are  highly  correlated.  Finally,  the  results  of  the  nonlinear  regression  from  all  77  profiles 
were  fit  to  a  multivariate  normal  distribution.  The  mean  and  standard  deviation  of  the  regression  parameters  are 
listed  in  Table  1.  Sample  realizations  of  these  fitted  profiles  are  plotted  in  Figure  6  along  with  the  prescribed 
profile.  The  measured  and  the  sampled  profiles  often  differ  significantly  from  the  prescribed  profile,  particularly 
near  the  edge  of  contact. 
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Figure  5.  Scatter  plots  of  the  profilometry  regression 
parameters,  y0  through  y12 


Figure  6.  Randomly  sampled  pad  profiles  compared  to  the 
prescribed  profile. 
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Table  1.  Summary  of  random  variable  statistics 


Random  V  ariable 

No. 

Mean 

St.  Dev. 

Distribution  Type 

Initial  crack  size,  a,  (in) 

A; 

5.95x1  O'4 

3. 34x1 0'4 

Lognormal 

Friction  coefficient,  juave 

a2 

0.302 

0.021 

Correlated  normal 

Partial  slip  slope,  a 

A, 

1.96 

0.120 

Pa  =  -  0.375 

Crack  growth,  log10(Cy) 

A, 

-12.7 

0.486 

Correlated  normal 

Crack  growth,  nt 

a5 

7.19 

0.715 

A  s  =  -0.9973 

Crack  growth,  log10(C2) 

a6 

-10.1 

0.157 

Correlated  normal 

Crack  growth,  n2 

x 7 

3.81 

0.146 

p67  =  -0.9751 

Profile,  y0  =  kL 

A« 

0.181 

5.84xl0'3 

Correlated  normal 

Profile,  yj  =  y0L 

a9 

-2335 

410 

Profile,  y2  =  RL 

A  w 

2333 

411 

Profile,  Y3  =  xoL 

A n 

-1612 

37.7 

Profile,  Yt  =  RR 

X]2 

2289 

379 

Profile,  Ys  = 

XI3 

1620 

35.4 

Profile,  Y6  =  kR 

XI4 

-0.183 

4.96xl0'3 

Profile,  y7  =  c7 

X 15 

-2.00xl0'4 

6. 20x1  O'4 

Profile,  yH  =  c2 

XI6 

-1.21xl0"6 

l.OlxlO'6 

Profile,  y,j  =  c3 

X 17 

1.53xl0"10 

6.38xl0'10 

Profile,  Yio  =  c4 

X 18 

9.80xl013 

6.16xl0'13 

Profile,  Yu  =  cs 

X 19 

-3.77xl0"17 

1.87xl0'16 

Profile,  Y12  =  c6 

x 20 

-3.80xl0"19 

1.59xl0"19 

IV.  Analysis  Methods 


A.  Monte  Carlo  Sampling 

Monte  Carlo  sampling  was  carried  out  to  determine  the  moments  (mean  and  standard  deviation)  of  the  cycles -to- 
failure.  This  procedure  involved  repeated  generation  of  realizations  of  the  random  variables  and  execution  of  the 
deterministic  fretting  fatigue  algorithm  to  determine  cycles-to-failure.  The  deterministic  fretting  fatigue  algorithm 
was  simplified  as  much  as  possible  to  minimize  the  computational  time  needed  for  each  set  of  random  variable 
samples.  This  included  calculation  of  the  contact  forces  using  the  method  demonstrated  by  Gean  and  Farris  [6], 
rather  than  modeling  the  contact  in  a  FEM  analysis.  The  CAPRI  run  time  was  minimized  by  reducing  the  number  of 
Fast  Fourier  Transform  terms  to  the  minimum  required  and  optimizing  the  solver  parameters  for  this  problem. 
Finally,  a  linear  fit  of  the  bulk  stresses  as  a  function  of  the  contact  forces  was  created  from  a  series  of  FEM  analyses, 
rather  than  using  a  new  FEM  analysis  for  each  Monte  Carlo  run.  Eliminating  the  FEM  analyses  from  the  calculation 
of  fretting  fatigue  life  reduced  the  computational  time  to  approximately  1.2  s  per  sample  when  running  in  parallel 
with  4  processors  on  an  Intel  Xeon  quad  core  workstation.  This  made  running  the  Monte  Carlo  analysis  with  a 
significant  number  of  samples  possible.  The  ensemble  of  cycles  to  failure,  Nf,  results  were  then  analyzed  to 
determine  the  mean  and  standard  deviation.  Sufficient  samples  were  executed  to  ensure  high  confidence  in  the 
computed  moments. 

B.  Sensitivity  analysis 

Probabilistic  sensitivities  play  an  important  role  in  determining  insight  into  the  dominant  factors  governing  a 
probabilistic  analysis  and  provide  information  as  to  potentially  effective  methods  to  improve  the  reliability.  There 
are  a  number  of  methods  in  the  literature  such  as  scatter  plots  [12],  Pearson  or  Spearman  correlation  [13],  regression 
methods  [14],  and  the  Score  function  method  [15],  among  others,  that  can  provide  useful  information. 

1 .  Scatter  plots 

Scatter  plots  are  a  two-dimensional  point  plot  of  the  sample  points  version  the  correspond  response  points.  The 
sample  realizations  for  each  random  variable  (X)  and  the  corresponding  results  (Y)  are  plotted  on  separate  axes.  If  a 
random  variable  is  not  important,  no  pattern  should  be  discernable,  that  is,  the  samples  for  X  should  mimic  its 
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marginal  distribution.  Conversely,  if  a  random  variable  is  important,  the  pattern  of  realizations  for  X  will  be 
distinctly  non-random.  Scatterplots  are  an  inexpensive  but  qualitative  method. 

2.  Regression  and  correlation 

Linear  regression  (LR)  methods  are  well  known  and  widely  available  tools  for  assessing  variance  contribution. 
LR  approximates  the  relationship  between  response  and  the  random  variables  as 

}’(x)  =  Ai  +  Z^/W  (6) 

where /is  an  arbitrary  function  of  the  random  variable  Xi.  In  general,  LR  may  contain  quadratic  and  interaction 
terms  of  any  number  of  variables.  The  correlation  coefficient  (Pearson  or  Spearman)  for  each  random  variable 
provides  an  indication  of  the  contribution  of  that  variable’s  variance  to  the  response  variance. 

3.  Segmented  input  distributions 

Sensitivities  based  upon  segmented  input  distributions  involve  dividing  the  input  samples  into  two  or  more 
separate  groups  based  on  the  output  samples.  The  output  and  related  input  samples  are  sorted  in  ascending  values 
based  upon  the  output  results.  If  the  variable  has  no  effect,  the  distributions  in  each  group  will  be  the  same  as  the 
marginal  distribution.  If  not,  then  the  variable  has  a  significant  effect.  Various  statistical  tests  such  as  Cramer-von 
Mises,  Mann- Whitney,  etc.  can  be  employed  to  determine  if  the  distributions  are  significantly  different. 

4.  Score  function 

The  Score  function  method  computes  the  partial  derivative  of  the  response  (cycles -to-failure)  with  respect  to  the 
parameters  of  the  input  PDFs.  The  probabilistic  sensitivities  (partial  derivatives  of  the  response  moments  (mean  and 
standard  deviation)  with  respect  to  the  parameters  of  the  random  variable  pdfs)  can  be  obtained  for  negligible  cost 
Probabilistic  sensitivities  such  as  0/tz/5/v,  represent  the  sensitivity  of  the  mean  of  cycles -to-failure  with  respect  to  the 
mean  of  random  variable  Xi. 


IV.  Results 


A.  Probabilistic  Life  Prediction 

The  probabilistic  fretting  fatigue  life 
prediction  tool  was  exercised  at  several 
experimental  load  cases.  Results  were  generated 
at  the  5  applied  loading  conditions  (Fmax  =18,  20, 
22,  24,  and  30  kN)  used  in  the  dovetail  testing 
program.  Lives  were  generated  at  each  loading 
condition  using  1000  samples  of  the  random  pad 
profiles,  coefficient  of  friction,  and  contact  force 
partial  slip  slope.  These  results  are  plotted  on 
lognormal  probability  paper  in  Figure  7  along 
with  the  experimental  lives.  The  experimental 
lives  shown  were  as  follows:  1,371,000, 
10,000,000,  and  693,000  cycles  at  18  kN; 
995,000  and  919,900  cycles  at  20  kN;  249,900, 
346,200,  and  497,800  cycles  at  22  kN;  164,000 
cycles  at  24  kN;  and  105,000  cycles  at  30  kN. 
These  data  were  quite  limited  for  an  evaluation  of 
a  probabilistic  life  prediction,  however,  they  were 
still  useful  to  ensure  the  predictions  were  near  the 
test  results.  The  comparison  in  Figure  7  shows 
that  at  the  higher  loads  the  predictions  were 
closer  than  at  the  lower  loads.  Since  the 
predictions  were  entirely  based  on  fatigue  crack 
growth  from  a  small  crack,  this  could  indicate 
that  a  significant  period  of  crack  formation  could 
be  taking  place  at  the  longer  lives. 

The  results  plotted  in  Figure  7  were  generated 
with  just  1000  samples  to  minimize  processing 


Figure  7.  Probability  plot  of  experimental  (points)  and 
predicted  (lines)  lives  at  several  applied  loads 
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time  for  each  condition  analyzed.  Since  the 
number  of  samples  was  relatively  small,  an 
understanding  of  the  variance  in  the  estimates  of 
mean  and  standard  deviation  of  Nf  was  important. 

Variance  estimates  of  the  modes  of  the  distribution 
were  calculated  at  each  of  the  applied  forces  and 
listed  in  Table  2.  Since  the  distribution  of  Nf  was 
fairly  lognormal,  the  modes  of  the  distribution  of 
the  logarithm  of  Nf  were  calculated.  The  columns 
showing  the  95%  confidence  bounds  of  the  mean 
were  then  converted  back  to  cycles.  The 
coefficient  of  variance  (COV)  was  calculated 
entirely  in  log  units.  Interestingly,  the  COV 
decreased  with  increased  applied  loading,  which 
was  consistent  with  fatigue  test  results  in  Ti-6A1- 
4 V  [10]  and  is  generally  true  in  metallic  aerospace 
materials.  Also,  in  Table  2  the  estimates  of  the 
modes  of  the  distribution  were  calculated  for  larger 
numbers  of  samples  at  Fmax  =  22  kN.  With  more 
samples,  the  confidence  bounds  shrink 
significantly  as  expected.  In  addition  to  quantifying  the  modes  of  the  distribution,  PDF’s  were  fit  to  the  distribution. 
Figure  8  is  a  plot  of  two  PDF  types  fit  to  the  Nf  results  at  Fmax  =  22  kN.  The  kernel  line  is  a  nonparametric  estimate 
of  the  PDF,  while  the  lognormal  curve  was  fit  using  the  mean  and  standard  deviation  of  the  logarithm  of  Nf.  The 


N  (cycles) 

Figure  8.  PDF  estimates  of  Nf  for  10,000  samples 


differences  in  the  PDF’s  showed  that  the  distribution  of  failure  lives  was  weighted  toward  longer  lives, 
also  be  observed  in  Figure  7,  and  this  trend  increased  with  lower  applied  loads. 
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Table  2.  Mean  and  standard  deviation  of  the  Monte  Carlo  analysis  results  including  confidence  bounds 

95%  confidence  for  mean  of  log  10(iV/)  (cycles)  95%  confidence  for  COV  of  log  10(iV/) 


Fmax  (kN)  Samples  Lower  Mean  Upper  Lower  COV  Upper 


18  kN 

1000 

833,200 

864,000 

895,900 

4.11% 

4.28% 

4.46% 

20  kN 

1000 

480,500 

495,600 

511,200 

3.66% 

3.82% 

3.98% 

22  kN 

1000 

319,400 

328,400 

337,700 

3.40% 

3.54% 

3.69% 

22  kN 

22  kN 

10,000 

50,000 

323,200 

326,200 

329,200 

3.67% 

3.72% 

3.77% 

24  kN 

1000 

215,600 

221,200 

227,000 

3.26% 

3.40% 

3.55% 

30  kN 

1000 

92,900 

95,100 

97,300 

3.08% 

3.21% 

3.35% 

B.  Sensitivity  Results 

1 .  Scatter  plots  and  correlations 

Figure  9  below  shows  the  scatter  plots  that  relate  cycles-to-failure  (y  axis)  with  each  random  variable  (x  axis)  for 
10,000  samples.  The  scatter  plot  for  random  variable  X2  shows  a  definite  pattern.  Correlation  coefficients  (Pearson 
or  Spearman)  for  each  variable  relate  to  the  amount  of  variance  in  Y  that  can  be  apportioned  to  X,.  The  results  for 
Pearson  correlation  coefficients,  given  in  Table  3,  indicate  that  variable  X2  (coefficient  of  friction)  is  by  far  the 
dominant  variable,  contributing  65%  of  the  variance  in  Nf.  It  must  be  pointed  out,  however,  that  the  correlation 
coefficients  are  obtained  without  consideration  for  the  simultaneous  variations  in  other  random  variables.  A  better 
estimate  of  variable  importance  can  be  determined  using  linear  regression,  discussed  below. 
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Figure  9.  Scatter  plots  of  the  model  output  cycles  to  failure,  Nfi  versus 
the  input  random  variables,  Xl  through  X20 


Table  3.  R-squared  values  for 
the  random  variables 


Variable 

R2 

x, 

0.024 

x2 

0.645 

x3 

0.093 

x4 

0.001 

X5 

0.002 

X6 

0.006 

Xy 

0.002 

X8 

0.002 

X9 

0.004 

Xjo 

0.004 

Xu 

0.003 

X]2 

0.005 

XI3 

0.017 

X,4 

0.008 

Xj5 

0.007 

X16 

0.063 

Xjy 

0.013 

Xis 

0.020 

X19 

0.020 

X2o 

0.038 

2.  Linear  Regression 

In  this  research,  an  LR  model  of  the  form 


y(X)  =  p0  +  YJPiXi 


(8) 


was  used,  where  X  denotes  a  vector  of  random  variables,  y  represents  the  cycles-to-failure  and  the  /?’ s  are 
coefficients  that  are  fit  to  the  analytical  results.  Table  4  shows  the  results  of  a  “best  model’’  fit  using  a  specified 
number  of  variables.  For  example,  using  only  2  variables,  the  combination  of  X2  and  X16  account  for  the  largest 
percentage  of  the  variance  in  Nf.  Note,  once  X2  is  included  in  the  model,  adding  X16  adds  17  percent  to  the  R2  sum, 
however,  the  R2  value  using  the  Pearson  correlation  coefficient  for  X2  without  considering  any  other  random 
variables  is  25%.  From  the  results  one  can  see  quickly  the  diminished  returns  offered  after  the  first  few  random 
variables.  For  example,  if  only  5  random  variables  are  used  in  an  LR  model,  this  model  would  account  for  85%  out 
of  a  possible  91%  of  the  output  variance.  It  is  somewhat  surprising  that  a  simple  linear  model  with  respect  to  the 
random  variables  can  account  for  such  a  large  percentage. 

Table  4.  Best  model  linear  regression  results 


- ^ - 

Number  R  C(p)  Random  Variables  in  Model 


1 

0.645 

28,370 

*2 

2 

0.713 

21,020 

*2 

X]6 

3 

0.754 

16,620 

V2 

X9 

Xio 

4 

0.824 

9000 

V2 

X]6 

Xis 

X2o 

5 

0.846 

6660 

X; 

X2 

Xi6 

Xis 

X2o 

6 

0.855 

5710 

xt 

X2 

X13 

XI6 

Xis 

X2o 

7 

0.874 

3640 

X ! 

X2 

x6 

Xy 

XI6 

Xis 

X20 

8 

0.883 

2690 

X , 

X2 

x6 

Xy 

Xu 

X,6 

X18 

X20 

9 

0.898 

1060 

Xj 

X2 

x6 

Xy 

X,6 

X2y 

X18 

Xj9 

10 

0.901 

720 

Xj 

X2 

x6 

Xy 

Xu 

XI6 

Xjy 

X18 

10 


Table  5  shows  the  results  using  LR  considering  natural  groupings  of  random  variables.  The  variables  were 
pardoned  into  4  groups:  (X,  -  initial  crack  size,  X2-X3  coefficient  of  friction  -  partial  slip  slope,  X4-X7  crack  growth 
parameters,  X8-X2o  geometry  profile).  The  results  indicate  that  group  2  (the  coefficient  of  friction-  partial  slip  slope 
is  dominant  followed  by  the  geometry  profile.  Surprisingly,  the  traditionally  dominant  random  variables  in  fatigue, 
crack  growth  rate  and  initial  crack  size  are  not  significant  in  terms  of  their  contribution  to  the  output  variance,  e.  g., 
these  parameters  could  be  modeled  deterministically. 


Table  5.  Linear  regression  group  analysis  showing  the  effect  of  each  random  variable 
group  on  the  model  R2 


Group 

RV  Added 

#  RV 

Model  R 2 

R1 

C(p ) 

F  value 

1 

X2,X3 

2 

0.645 

0.645 

28,370 

9080 

2 

Xg  —  X2  o 

15 

0.210 

0.855 

5710 

1110 

3 

X4-Xy 

19 

0.030 

0.885 

2420 

660 

4 

X , 

20 

0.022 

0.908 

20 

2410 

3.  Segmented  Distributions 

10,000  samples  were  sorted  by  cycles-to-failure  then  segmented  into  three  groups:  early  failures,  middle 
failures,  and  late  failures.  That  is,  samples  in  the  first  group  failed  early  (low  cycles-to-failure),  and  samples  in  the 
third  group  failed  last,  etc.  Nonparametric  density  estimation  methods  [16]  were  used  to  develop  PDF’s  for  each 
random  variable  for  each  group  of  samples  and  for  the  total  ensemble  of  samples.  Figure  10  shows  the  PDF’s  for 
three  random  variables  that  exhibit  a  large,  moderate  and  minimal  effect  due  to  segmentation.  The  violet,  red,  and 
green  distributions  derive  from  the  early,  middle  and  late  cycles-to-failure.  The  black  dashed  line  indicates  the  PDF 
using  the  total  ensemble  of  samples.) 


Figure  10.  Segmented  PDF’s  showing  large  (X2),  moderate  (XI6)  and  minimal  (Z5)  effects 


V.  Summary  and  Conclusions 

A  new  probabilistic  fretting  analysis  was  developed  to  investigate  the  relative  importance  of  typical  fretting  input 
variables  on  the  predicted  failure  lives.  Several  random  variable  inputs  were  identified  and  pdfs  were  quantified 
using  laboratory  data.  Monte  Carlo  sampling  of  the  input  pdfs  was  performed  and  a  deterministic  analysis  was 
repeatedly  run  using  the  sampled  inputs  to  obtain  a  distribution  of  predicted  fretting  lives.  The  results  showed  that 
significant  scatter  in  fretting  lives  can  be  expected  based  on  variability  in  the  material  properties,  contact  profiles, 
coefficient  of  friction,  and  contact  force  response.  Interestingly,  the  dominant  variables  in  terms  of  contribution  to 
the  cycles-to-failure  variance  were  the  coefficient  of  friction  and  several  terms  within  the  geometry  profile;  whereas, 
the  traditionally  dominant  variables,  initial  crack  size  and  crack  growth,  were  not  significant. 
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